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ABSTRACT 

•  ■  Vcrj; 

We- propose  a  hybrid  method  for  computing  the  feedback  gains  in  linear 
quadratic  regulator  /problems.  The  method,  which  combines  use  of  a  Chandrasekhar 
type  system  with  an  iteration  of  the  Newton-Kleinman  form  with  variable 
acceleration  parameter  Smith  schemes,  is  formulated  so  as  to  efficiently  compute 
directly  the  feedback  gains  rather  than  solutions  of  an  associated  Riccati  equation. 
The  hybrid  method  is  particularly  appropriate  when  used  with  large  dimensional 
systems  such  as  those  arising  in  approximating  infinite  dimensional  (distributed 
parameter)  control  systems  (e.g.,  those  governed  by  delay-differential  and  partial 
differential  equations).  Computational  advantages  of  our  proposed  algorithm  over  the 
standard  eigenvector  (Potter,  Laub-Schur)  based  techniques  are  discussed  and 
numerical  evidence  of  the  efficacy  of  our  ideas  presented. 


Kev  Words:  LQR  problems,  feedback  gains,  distributed  parameter  systems, 
computational  algorithm,  Chandrasekhar  system,  Newton-Kleinman  scheme,  Smith 
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A  great  deal  of  effort  in  recent  years  in  control  of  distributed  systems  has 


focused  on  approximation  techniques  (for  example,  see  [l]-[6],  [8],  [1 1  ]-[  1 4],  [16],  [18], 


[21],  [22],  [24],  [26],  [32])  to  reduce  inherently  infinite  dimensional  problems  to  (large) 


finite  dimensional  analogues.  Relatively  little  attention  has  been  given  to  the 


development  of  efficient  computational  methods  for  the  resulting  large  but  finite 


dimensional  control  problems.  In  this  paper  we  consider  such  questions  for  one 


classical  formulation  of  the  feedback  control  problem,  the  well-known  linear  quadratic 


regulator  (LQR)  problem. 


There  are  several  approaches  one  can  take  in  such  an  endeavor.  With  the 


emergence  of  new  computer  architectures  (vector  and  parallel),  one  exciting  possibility 


involves  the  development  of  new  algorithms  to  be  used  with  nonsequential  computers. 


While  we  are  currently  investigating  ideas  in  this  direction,  our  presentation  here  is 


to  report  on  some  of  our  efforts  to  develop  better  algorithms  for  use  with 


conventional  serial  computers. 


As  is  well-known,  the  LQR  problem  can  be  reduced  to  the  solution  of  a 


matrix  Riccati  equation  in  order  to  construct  the  feedback  gain  matrix.  The  most 


widely  available  method  for  solution  of  the  Riccati  equation  is  the  Potter  method 


[30],  which  is  based  on  obtaining  the  eigenvectors  and  eigenvalues  of  an  associated 


2nx2n  Hamiltonian  system  when  the  underlying  dynamical  control  system  is  of 


dimension  n.  A  related,  but  improved  version  involving  computation  of  Schur  vectors 


for  the  system  was  proposed  by  Laub  in  [27],  While  both  of  these  "eigenvector" 


methods  can  be  used  satisfactorily  (for  a  discussion  of  real  advantages  offered  by  the 


Laub-Schur  approach  over  Potter’s  method,  see  [27])  for  systems  with  n  relatively 


small,  say  n  <  100,  the  computational  effort  (and  time)  grows  like  n*  and  becomes 


prohibitive  for  large  systems.  More  recently,  the  idea  of  using  Chandrasekhar  systems 


((9],  [20],  [33,  p.304-310],  [36])  when  the  number  of  states  is  large  compared  to  the 
number  of  control  inputs  (exactly  the  situation  in  a  number  of  cases  where  one 
approximates  a  distributed  system)  has  been  suggested  by  a  number  of  authors  [33, 
p.309],  [7],  [8],  [17],  [31].  However,  as  we  shall  discuss  below,  there  can  be  numerical 
difficulties  in  using  the  Chandrasekhar  approach.  On  the  other  hand,  it  is  known 
that  iterative  methods  such  as  the  Newton  algorithm  as  formulated  by  Kleinman  in 
[23]  can  be  quite  efficiently  implemented  (even  for  some  large  systems)  if  good 
initial  estimates  are  provided  and  if  one  can  solve  efficiently  the  resulting  Lyapunov 
equations.  In  this  presentation,  we  discuss  the  formulation  and  numerical  testing  of  a 
hybrid  method  that  represents  an  attempt  to  combine  the  good  features  of  the 
Chandrasekhar  approach  (growth  like  n  in  computational  effort)  with  those  of  the 
Newton-Kleinman  (quadratic  convergence  when  good  initial  estimates  are  provided) 
along  with  innovative  use  of  the  Smith  algorithm  for  solution  of  Lyapunov 
equations. 

We  expect  these  ideas  to  be  quite  useful  in  design  of  control  laws  for  some 
of  the  models  currently  being  investigated  in  connection  with  large  flexible  structures 
as  well  as  in  some  of  the  population  dispersal  and  control  studies  that  we  are 
currently  pursuing  with  biologists  and  ecologists.  Some  of  the  large  flexible  structures 
involve  rather  sophisticated  distributed  parameter  models  (e.g.  see  [4],  [34]),  especially 
when  one  wishes  to  include  complicated  damping  mechanisms  involving  time  or 
spatially  related  hysteresis  [6],  [34]  or  nonlinear  effects  [19].  For  such  models,  the 
computational  task,  can  be  rather  demanding  whether  one  is  carrying  out  parameter 
identification  [3]  or  feedback  control  calculations  with  traditional  eigenvector  based 
methods  (the  authors  of  [6]  have  indicated  experiences  with  runs  requiring  9  hours  of 
VAX  time  when  using  an  approximate  system  with  dimension  equivalent  to  n  « 
238). 

For  our  presentation,  we  assume  that  one  has  used  their  favorite 
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approximation  scheme  (finite-elements,  spline,  spectral,  etc.)  to  reduce  the  problems  of 
interest  to  an  LQR  problem  with  finite  dimensional  system.  More  precisely, 
throughout  our  discussions  we  consider  the  LQR  problem:  minimize  the  cost 
functional 

(1.1)  J(u)  -  [*{  |Cx(t)|2  +  |u(t)f}dt 

^  o 

subject  to  the  state  dynamics 

x(t)  -  Ax(t)  +  Bu(t)  ,  x(0)  -  x0  . 

Here  A  €  RnXn,  B  €  RnXm,  and  C  e  RpXn.  (We  have,  without  loss  of  generality,  for 
our  discussions  here  normalized  our  problem  so  that  the  control  term  in  the  cost 
functional  (1.1)  appears  with  a  weighting  matrix  I.)  We  shall  assume  that  (A,B)  is 
stabilizable  and  (C, A)  is  detectable  [25],  [37].  Then  the  optimal  feedback  control  for 
the  LQR  problem  involving  (1.1)  is  given  by 
u(t)  =  -BTPx(t) 

where  P  is  the  unique  non-negative  symmetric  solution  of  the  algebraic  Riccati 
equation 

(1.2)  AtP  +  PA  -  PBBtP  +  CTC  -  0. 

In  this  paper  we  propose  an  algorithm  which  leads  to  direct  calculation  of  the 
feedback  gain  matrix  K  =  BTP  without  computation  of  P.  In  addition  to  providing 
substantial  savings  in  computational  time  over  eigenvector  methods,  our  algorithm 
requires  much  less  storage  and  can  easily  be  implemented  to  take  advantage  of 
special  structures  (e.g.  sparsity)  in  the  system  matrix  A. 

To  outline  the  steps  in  our  algorithms,  we  first  recall  that  the  optimal 
feedback  gain  K  is  given  by  the  limit  K  ■  lim  K(t)  as  t  -  •*  of  solutions  of  the 
Chandrasekhar  system  [9],  [20] 


(1.3) 


K(t)  -  -BTLT(t)L(t)  ,  K(0)  -  0 
L(t)  -  -L(t)(A-BK(t))  ,  L(0)  «  C, 
where  K  €  RmXn  and  L  e  RpXn.  In  fact  (see  [20],  [37])  P  =  lim  J°  LT(s)L(s)ds  as  t  - 
-*».  The  first  step  in  the  proposed  hybrid  algorithm  involves  a  numerical  integration 
of  (1.3)  backward  in  time  on  an  appropriate  interval  [-t^O].  For  the  second  step,  the 
value  K(-tf)  obtained  through  this  numerical  integration  of  the  Chandrasekhar  system 
is  then  refined  by  use  of  the  Newton-Kleinman  algorithm  [23],  if  we  use  K(-tf)  as  an 
initial  value  K0  for  the  Newton  method. 

To  motivate  our  effort  in  these  two  steps,  we  note  that  the  convergence 

K(t)  "  K  as  t  -  •*  can  be  very  slow  when  the  eigenvalues  of  A-BK  lie  close  to  the 

imaginary  axis.  Moreover,  L=0,  K*  arbitrary  are  solutions  in  the  asymptotic  limit 
sense  to  (1.3).  That  is,  if  we  denote  by  f(K,L)  the  right  side  of  system  (1.3),  then  Km 
arbitrary  and  L=0  are  solutions  of  f(K«,L)  =  0.  Hence  K(t)  •*  K»,  L(t)  -  0  as  t  -*  -» 
doesn’t,  in  general,  have  a  unique  limit  numerically.  Thus,  as  is  pointed  out  in  [33,p. 
316-318],  if  one  is  to  use  the  Chandrasekhar  approach  alone,  one  needs  a  very 
accurate  numerical  solver  for  (1.3).  This  can  be  computationally  quite  expensive  if  we 
are  dealing  with  a  large  system  and/or  a  stiff  system.  Hence,  we  propose  to  use  a 
rather  crude,  fast  integration  method  for  the  Chandrasekhar  component  of  our 
algorithm  and  take  the  resulting  numerical  solution  K(-tf)  as  a  start-up  value  for  the 
Newton  iterations.  If  this  crude  estimate  from  the  Chandrasekhar  step  is  a 

sufficiently  good  initial  guess,  then  we  can  expect  to  meet  the  Newton-Kleinman 
requirements  that  A-BK0  be  a  stability  matrix  and  to  obtain  quadratic  convergence  in 
this  second  component  of  the  algorithm. 

The  first  step  of  our  hybrid  method  requires  the  solution  of  n(m+p) 

simultaneous  equations,  while  each  iteration  of  the  usual  Newton-Kleinman  step 
requires  the  solution  of  a  Lyapunov  equation  for  the  nxn  symmetric  estimates  of  P. 


However,  as  we  shall  see  below,  one  can  use  factorization  ideas  [20]  and  the  Smith 
method  [35]  for  Lyapunov  equations  to  reformulate  the  Newton-Kleinman  method  as 
a  direct  iterative  method  for  the  mxn  gain  K,  thereby  providing  additional 
computational  advantages.  To  speed  up  our  calculations  and  improve  convergence  in 
the  Smith  algorithm,  we  propose  a  variable  stepsize  Smith  method  to  solve  the 

Lyapunov  equations  as  described  in  Section  4  below.  In  Section  2  we  outline  a 

numerical  scheme  for  the  Chandrasekhar  system,  while  the  reformulated 

Newton-Kleinman  iterative  procedure  to  compute  directly  the  gain  K  is  detailed  in 

Section  3.  Finally,  in  Section  5,  we  discuss  further  some  advantages  and  disadvantages 
of  the  proposed  algorithm  and  report  on  our  experience  with  several  numerical 
examples  to  illustrate  the  feasibility  of  our  hybrid  approach. 


2.  Numerical  Solution  of  the  Chandrasekhar  System 


Wc  return  to  consider  more  closely  the  Chandrasekhar  system: 

(2.1)  K(t)  -  -BTLT(t)L(t)  ,  K(0)  -  0 

(2.2)  L(t)  *  *L(t)(A  -  BK(t))  ,  L(0)  -  C  . 

where  A  t  Rn*n,  L,C  t  Rp*n,  and  K,BT  t  RmXn.  We  first  observe  that  the  second 
equation  (2.2)  is  linear  in  L.  In  cases  where  A  arises  from  a  discretization  or 
approximation  of  partial  differential  equations,  equation  (2.2)  tends  to  be  a  stiff 
system  and  thus  it  is  advisable  to  use  an  implicit  numerical  scheme.  We  propose  here 
the  second  order  Adams-Moulton  algorithm  [10,  p.235],  A  second  observation  is  that 
the  right  hand  side  of  equation  (2.1)  is  independent  of  K  and  thus  an  explicit 
scheme  is  appropriate;  we  propose  the  second  order  Adams-Bashforth  algorithm  [10, 

p.226]. 

These  observations  lead  us  to  propose  the  following  algorithm  for  the 
Chandrasekhar  system  (2.1),  (2.2):  Given  a  step  size  h>0,  approximations  K;  and  L;  to 
K(-ih)  and  L(-ih)  are  generated  by 


(2.3) 

-  Ki  + 

h(3-BTLT  Li  - 

1-BTLT  L 

(2.4) 

•Cs 

(0) 

-  <1 

+  Kj)/2 

(2.5) 

^i+i 

**  Li  + 

j<Li«  +  MHa-bk'^) 

(2.6) 

"  Kj  4 

where  K0  - 

0  and  L_j 

=  Lo  - 

c. 

Several  remarks  may  be  useful  at  this  point. 

(Remark  1)  The  stiffness  of  the  matrix  A  dictates  the  choice  of  stepsize  h. 
(Remark  21  The  predicted  values  |  and  the  corrected  values  Ki+1  satisfy 
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Ki+1  '  Ki+1  “  |BT(LH-lLi+X  ‘  2L^Li  +  LUlLi-l) 

“  I  BT  £ 

and  this  relationship  can  be  used  for  stepsize  refinement,  i.e.,  to  give  local  bounds 
depending  on  stepsize  which  can  be  used  in  error  control. 

(Remark  3)  The  formula  (2.5)  can  be  rewritten  as 


(2.7)  Li+1  =  Lj(I  +  |  A|)(I  -  hA.)-t 

=  2Lj(I  -  jAj)1  -  Lj 

where  A.  s  A  -  Bk|°^  .  Defining  H  ■  I  -  jlA  we  have  that  I  -  ^Aj  ■  H  + 
where  B  €  RnXm  and  e  RmXn  .  Thus  by  the  Sherman-Morrison-Woodbury 

formula  [29,  p.50]  (used  frequently  when  updating  an  nxn  matrix  by  rank  m 
matrices) 


(2.8) 


(I  -  ^A;)'1  -  H-1  - 


where  I  +  €  Rm*m.  The  matrix  K^^H*1  in  (2.8)  can  be  computed  by 


(2.9) 


*!>-'  -  KH-1  *  ^|BTL*Lp-‘  - 
K1+1H-*.  K,H->  +  ^BTLj+,L1+1H->  ♦  B'LfL.H-'l. 

Hence  we  see  from  (2.7)  -  (2.9)  that  the  step  (2.5)  only  involves  the  operation  LjH'1 
plus  inversion  of  an  mxm  matrix  I  +  H^B.  Thus  the  step  (2.5)  can  be 

reformulated  so  that  it  requires  only  an  mxm  matrix  inversion  plus  matrix-vector 


multiplications  if  the  LU  decomposition  of  H  *  I  •  |a  is  computed  a  priori.  This 


procedure  can  be  most  advantageous  computationally  when  m  and  p  are  small 
compared  to  n. 

(Remark  4)  For  some  problems  one  might  wish  to  use  a  completely  implicit  scheme 
in  place  of  (2.3)  -  (2.6)  to  enhance  stability  and  reduce  sensitivity  to  step  size  choice. 
Then  one  might  consider  iterations  Kp+[  ,  l|j+{  generated  by 


£ 


(2.10) 
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Li+i  "  Li  +  ^l!+i  +  Li)(A-BK [$) 

K.+i  -  Ki  +  ♦  bTlTls) 

„.0)  ,  (j) 

"  (^i+i  +  ^i)/2 


and  thus  produce  iterates  with  limits  (as  j  -  •)  K;+1,  Ll+1  satisfying 


(2.11) 


Li+i  “  Li  +  5<Li+i  +  Li)(A-B(Ki+1  +  Kj)/2) 


K: 


T  T 


T  T 


i+1 


K,  ♦  S(B  Li+IL1+1  ♦  B  L,  L,) 


r>. 


A  widely  used  iterative  method  for  finding  the  non-negative  solution  of  the 
algebraic  Riccati  equation  (1.2)  is  Newton’s  method  as  modified  by  Kleinman  [23].  We 
show  that  this  method  can  be  reformulated  so  that,  when  combined  with  a  factored 
form  of  the  well  known  Smith  method  [35],  one  can  compute  directly  a  sequence  of 
iterates  Kj  for  the  feedback  gain  K. 

First  we  recall  the  Newton  iterative  algorithm  as  formulated  by  Kleinman: 

(1)  Choose  a  gain  matrix  K0  so  that  A-BK0  is  a  stability  matrix; 

(2)  Update  K;  by  Ki+1  =  BtPj  where  P;  is  the  solution  of  the  Lyapunov  equation 

(A-BK^Pj  +  Pj(A-BKj)  +  kTK;  +  CTC  -  0. 

It  can  be  shown  that  0  S  Pi+1  <  Ps  for  any  i,  and  P  =  lim  Pj  where  the 
convergence  is  quadratic.  This  algorithm  can  be  viewed  as  an  iterative  method  for 
the  gain  K,  i.e.  K  =  lim  Ks  where  Ki+1  »=  F(K;)  with  F(K)  =  BTP  and  P  is  the 
solution  of  the  Lyapunov  equation: 

(3.1)  (A-BK)tP  +  P(A-BK)  +  KtK  +  CTC  -  0. 

Thus,  in  order  to  calculate  F(K)  one  has  to  solve  the  Lyapunov  equation  (3.1)  for  the 
symmetric  matrix  P.  However,  one  can  form  an  alternative  version  that  allows  one  to 
directly  calculate  F(K)  using  the  Smith  method  for  a  Lyapunov  equation  in  X  of  the 
form 

(3.2)  STX  +  XS  +  DtD  =  0 

where  S  €  RnXn  i$  a  stability  matrix  and  D  €  RpXn. 

To  this  end,  we  replace  step  (2)  in  the  Newton-Kleinman  method  by: 
(2')  For  i  >  1,  update  Kj  by  Ki+1  =  K|  ■  BtZj  where  Z,  ■  PU1  -  Pj  is  the  solution  of 
the  Lyapunov  equation 


-10- 


with  Dj  s  Kj  •  Kj_r 

The  method  with  (2')  offers  several  advantages  over  that  using  (2).  The 
Lyapunov  equation  in  (21)  has  fewer  inhomogeneous  terms  than  does  the  one  in  (2) 
and  the  term  D;  has  rank  m  which  depends  only  on  the  number  of  inputs  (controls) 
to  the  system.  In  the  proposed  modified  Smith  method  described  below,  one  is  able 
to  compute  directly  the  m*n  update  matrix  J1  =  BTZ;  without  computing  Z.  (see  (3.12) 
below).  Since  Z;  -*  0  as  i  -*  •  is  expected,  choosing  the  start-up  value  JJ,  «  0  in  the 
factored  Smith  algorithm  (where  J1  ■  BTZj  is  computed  as  the  limit  as  k  -*  •  of  a 

sequence  Jj^)  is  a  natural  as  well  as  convenient  choice. 

Note  that  the  step  (2*)  requires  that  one  have  Zj  *  P0  -  P2  in  hand  and 
hence  we  must  start  this  procedure  with  P0,  Pt  (and  K0,Kj)  given  whereas  (2)  requires 

only  that  one  start  with  K0  given.  Then  Kx  is  computed  by  Kx  -  BTP0  with  P0  the 

solution  of 

(A-BK0)tP0  +  P0(A-  BK0)  +  Kj  K0  +  CTC  =  0. 

Since  our  Smith  algorithm  below  is  formulated  to  solve  Lyapunov  equations  of  the 
form  (3.2),  we  can,  to  maintain  this  form,  initially  solve  the  equation  twice.  That  is, 
if  we  solve  for  Z0  the  solution  of 

(3.4)  (A-BK0)TZ0  +  Z0(A-BK0)  +  K*  K0  «  0 
and  for  Z0  the  solution  of 

(3.5)  (A-BKq)TZ0  +  Z0(A-BK0)  +  C* C  «  0  , 

then  we  can  obtain  Kj  by  =  BTZ0  +  BTZ0  .  Since  the  Smith  method  as 
formulated  here  actually  returns  BTX  where  X  is  the  solution  to  (3.2),  we  thus  will 
use  this  Smith  algorithm  twice  (with  S  *  A-BK0),  once  with  DTD  ■  KjK0,  once  with 
DtD  -  CTC  and  then  simply  add  the  solutions  to  obtain  Kr 

We  turn  next  to  the  desired  factored  form  of  the  Smith  method  as  applied 
to  equation  (3.2).  Let  X0  be  an  arbitrary  n*n  symmetric  matrix  and  let  a  sequence 


{Xj,}  of  nxn  symmetric  matrices  be  generated  by 

(3.6)  Xk+1  -  UjX.U,  +  Yr  , 

where  r  is  a  postive  constant  (the  Smith  stepsize)  and 

(3.7)  Ur  -  (I-rS)'1  (I+rS) 

(3.8)  Yr  -  2r(I-rST)"1  DTD(I-rS)*1  . 

Then  one  can  argue  that  (Xk)  converges  to  X,  the  solution  of  (3.2).  The  method  and 
its  analysis  is  based  on  the  observation  that  for  any  positive  constant  r,  the  equation 
(3.2)  is  equivalent  to 

X  «  U*  XUr  +  Yr 

which  can  be  used  to  define  a  contraction  map  in  the  obvious  manner. 

We  modify  this  standard  formulation  of  the  Smith  method  to  suit  our 
particular  needs  here  (computing  BTX  instead  of  X).  From  (3.6)  we  have 

-  Xk  -  U?<Xk  •  Xk.,)U,  .  k  >  1. 

Hence,  if  Xk  -  Xkl  =  MkMk  (i.e.  if  we  have  a  factorable  difference),  then 

Xk+x  -  Xk  “  U7  Mk  MkUr  “  (MkUr)T(MkUr). 

If  the  start-up  value  X0  is  zero,  then  we  can  write 

(3.9)  Xt  -  X0  =  2rMj‘M1  ,  Mx  *  D(I-rS)’1  . 

By  induction  on  k,  (3.6)  is  then  equivalent  to 

(3.10)  Mk+1  =  MkUr 

(3.11)  Xk+1  «  Xk  +  2rMk+1Mk+1  . 

In  this  manner  BTX  can  be  obtained  as  the  limit  of  Jk  -  BTXk  where  Jk  is  generated 
by 

Jk+1  -  Jk  +  2rBX+iMk+1  • 

Thus,  the  update  step  (21)  is  carried  out  by  the  following  Smith  algorithm: 

(3.12) 0)  Set  Sj  -  A-BKj  and  D  ■  K.  -  KU1; 

(3.12) (»)  Choose  a  positive  constant  r  and  form  Ur  and  Mj  by  (3.7)  and  (3.9) 
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with  S  -  Sj  ;  put  J0  *  0  and  Jj  -  J0  +  2rBTM^Mr 
(3.12)(m)  Iterate  for  k  =  1,  2,  .  .  .  ,  on 

Mk+1  =  MkUr 

Jktl  -  Jk  +  2rBTMk+1Mk+1  . 

In  summary,  we  have  described  in  this  section  a  Newton-Kleinman  scheme 
combined  with  the  Smith  method  for  the  resulting  Lyapunov  equation  at  each  step  in 
the  Newton-Kleinman.  We  have  reformulated  the  Newton-Kleinman  iteration  and 
factored  the  Smith  algorithm  so  as  to  result  in  algebraic  savings  in  computing 
directly  the  gain  estimates  Kj. 


As  is  well-known,  the  rate  of  convergence  in  the  Smith  method  discussed  in 
the  last  section  depends  upon  the  choice  of  the  acceleration  or  step  parameter  (See 
[33,  p.291-297]  for  several  discussions.  Note  that  our  parameter  r  is  the  negative 
reciprocal  of  the  parameter  in  Russell’s  discussions).  To  increase  speed  in 
convergence,  one  may  employ  the  accelerated  Smith  method  [33],  [35]  which  can  yield 
quadratic  convergence  as  compared  to  the  linear  convergence  obtained  with  (3.6). 
However,  unlike  (3.6),  the  accelerated  Smith  method  is  not  self-correcting  [33]  and 
here  we  propose  to  speed  up  convergence  in  an  alternative  way  which  has  proved 
both  reliable  and  efficient  in  some  of  our  numerical  tests.  Specifically,  we  propose  to 


use  a  succession  of  acceleration  parameter  values  rt  (much  in  the  spirit  of  other 
well-known  iterative  methods  such  as  alternating  directions  [15],  [28])  to  accelerate 
convergence  in  the  basic  Smith  method.  Our  formulation  of  this  "variable  stepsize" 
Smith  method  is  based  upon  the  observation  that  for  fixed  r  >  0  and  k  >  1,  the 
Smith  algorithm  can  be  written  as 

STxk  +  Xks  +  DtD  -  Ek 

(41)  T  T 

Ek  s  (I  +  rS  )TMk  Mk(I  +  rS) 

where  Mk  is  defined  by  (3.9)  and  (3.10).  To  see  this,  we  note  that  from  (3.6)  we 


(I  -  rS)TXk(I  -  rS)  -  (I  +  rS)TXk  l(I  +  rS)  +  2rDTD  , 


(I  +  rST)Xk(I  +  rS)  -  (I  -  rS)TXk(I  -  rS)  +  2rDTD 


(I  +  rS)T(Xk  -  Xk.jKI  +  rS)  . 


Hence,  from  (3.11)  we  obtain 


2r(STXk  +  XkS  +  DtD)  -  2r(I  +  rS)TMjMk(I  +  rS) 


which  implies  (4.1).  Moreover,  from  (3.9)  and  (3.10)  we  have 
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(4.2)  Ek  -  (Uj  )kDTD(Ur)k  ,  k  >  1  . 

Thus,  if  we  use  the  iteration  (3.6)  with  acceleration  parameter  rx  for  kj  iterates,  we 
obtain  an  iterate  X'1^  and  equation  error  satisfying 


STX(1)  +  X(%  +  DtD  -  E(1) 


,  ,  k  k 

E(1)  -  (U?)  1  DtD(U  )  l. 

rt  ri 

Let  us  define  the  difference  E^  -  X  -  X^  where  X  is  the  sought-after  solution  of 


(3.2).  Then  it  is  readily  seen  that  E^  satisfies  a  Lyapunov  equation  similar  to  that 


of  (3.2)  : 


STI  +  IS  +  E(1)  «  0. 


If  we  next  apply  the  iteration  (3.6)  k,  times  with  acceleration  parameter  r,  to  the 


residual  equation  (4.4)  we  obtain 


STX(2)  +  X(2)S  +  E(1)  -  E{2> 


where  X^  is  the  final  iterate  using  r2  and  the  equation  error  E^  is  given  by 


E(2)  -  (Uj)kj  E(1)(Ur  )kj  . 
r2  r2 

If  we  proceed  to  define  the  difference  E*2)  -  X  -  (X*1*  +  X(2^),  then  from  (4.4)  and 
(4.5)  we  sec  that  E^  satisfies  a  Lyapunov  equation 


StE  +  ES  +  E(2)  -  0  . 


We  continue  this  procedure,  using  a  sequence  of  acceleration  values  {rs)  along 
with  corresponding  iteration  counts  {kj}  to  produce  a  sequence  {X^}  of  nonnegativc. 


symmetric  matrices.  For  i  )  1,  we  have 


STX«  +  X(i)S  +  E(i'1}  -  E(i) 


(4.7)  E(i)  -  (U£)ki  E(i-1J(Ur)ki  ,  E*°)  e  DTD  . 

Thus,  if  Xj  *  £  x(i),  then  Xj  satisfies 


STXj  +  XjS  +  DtD  -  E(j) 


and  hence  X;  <  X  and  X,  ,  <  X,  ,  j  >  1. 
j  j 


Using  arguments  similar  to  those  in  [33,  p.291-297]  one  can  show  that  for 


0  <  £  <  r  <  R,  with  £,  R  positive  constants,  there  exists  a  constant  u,  0  <  u  <  1, 
depending  only  on  £  and  R,  such  that  for  u  <  p  <  1, 


|Uk  |  <  M(p)pk  ,  k  -  0,  1 . 

where  M(p)  is  independent  of  r.  Thus,  if  £  <  r{  <  R,  then  for  any  0  <  6  <  1,  there 
exists  an  integer  k(6)  such  that  for  kj  >  k(5)  ,  i  >  1, 

(4.9)  |U*‘|  <1-6. 

Hence,  using  (4.7)  we  have 

|  E<j)  |  <  (l-£)2j  j  D  J 3 

so  that  E^  -«  0  as  j  -  •  and  therefore  X,  •*  X  as  j 

For  the  hybrid  method  proposed  in  this  paper,  we  have  combined  the 
variable  stepsize  method  just  outlined  with  the  reformulated  Smith  method  of  (3.12). 
We  then  obtain  the  following  algorithm  for  solving  for  the  feedback  gains  Kt: 


Algorithm  (4, IQ): 

Set  Sj  -  A  -  BK.,  Dj  -  Kj  -  KU1  and  J0  -  0.  For  given  acceleration  parameters  rv 
r2,  .  .  .  ,  and  iteration  count  indices  kp  k2,  ....  we  iterate  on  j  *  1,  2,  ....  in 
the  following  steps: 

(4.10a)  Compute  Ur  -  (I  -  rjSi)_1(I  +  r^S,)  and 

“i  *  Dfl  -  Wl  . 

-  J0  +  2rjBTM*MI 

(4.10b)  Iterate  for  k  -  1,  2,  .  .  .  ,  k  j  -  1  in 

Mk+1  «  MkUr. 

V,  -  J„  -  2'iBTMf+iMk« 

(4.10c)  Compute  Dj+1  -  Mk  (I  +  r^Sj),  set  J0  ■  Jk  and  return  to  (4.10a)  with 


M 

N 
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The  steps  in  Algorithm  (4.10)  are  to  be  repeated,  i.e.  iteration  through 
r.,  r 2,  .  .  .  ,  until  some  convergence  criterion  is  met.  In  performing  the  steps  in 
(4.10b),  one  can  use  the  Sherman-Morrison-Woodbury  formula  in  a  procedure  such  as 
that  outlined  in  (2.7)  -  (2.9)  in  section  2. 

As  is  often  the  case  in  variable  stepsize  algorithms,  the  choice  of  the 


acceleration  parameters  rr  r2,  .  .  .  ,  and  the  associated  iteration  counts  kj,  k2. 


provides  both  freedom  and  some  frustration  in  the  search  for  "best"  choices.  If  one 
follows  the  guide  provided  by  ADI  methods  (see  [15,  p.37]),  one  might  choose  a  set 
of  values  r.  to  be  used  in  some  cyclic  order.  The  best  choices  of  values  for  the 
often  depend  on  the  eigenvalues  of  Sj  *  A  -  BK;.  For  example,  consider  the  case 
where  S4  has  only  real  eigenvalues  X^,  each  with  multiplicity  mj,  j  «  1,  2,  ....  m. 
Then  a  choice  of  r.  =  -I/X.  and  kj  =  m^  in  the  algorithm  produces  convergence  in  a 
finite  number  (m)  of  steps.  That  is,  this  choice  yields  «  0  in  (4.8). 

Of  course,  the  complete  eigenstructure  of  S;  will  not  be  known  (nor  do  we 

suggest  that  any  sophisticated  analysis  along  these  lines  be  included  with  each  use  of 

Algorithm  (4.10)  to  obtain  the  gains  K;).  A  possible  alternative  is  to  use  one  of  the 

polynomial  acceleration  methods  [15,  Chp  3,  4], 

In  closing  this  section,  we  note  that  the  analogies  of  our  variable  step  Smith 
method  with  the  ADI  methods  used  to  solve  partial  differential  equations  can  be 
made  a  little  more  precise.  Briefly,  in  ADI  spliting  methods  [28,  p.146-148],  one 
attempts  to  solve  a  discretization  of  the  evolution  equation 

4>  *  A$>  +  f 

when  A  >  0  can  be  written  A  =  Aj  +  A2  with  Aj  >  0  (for  example,  factored  into 
components  corresponding  to  spatial  discretizations  in  the  x  and  y  directions 
respectively  for  an  equation  in  a  two  dimensional  spatial  domain).  This  can  be  shown 
[28,  p.  1 50]  to  lead  to  an  iterative  scheme 


(4.11) 


(1  +  |Aj  )(i  +  ^a^*1  -  (i  -  ^AjKi  -  |a2)<^  +  hfj 


where  the  index  j  is  related  to  time  stepping.  On  the  other  hand,  if  one  considers 
the  Smith  method  (3.6)-(3.8)  for 
STX  +  XS  -  F 


and  chooses  r  =  li,  one  obtains  the  iteration 
2 


(4.12) 


(I  -  ^ST)Xj+1(I  -  |S)  -  (I  +  £sT)Xj(I  +  tS)  +  hF  . 


In  these  iterations  one  may  identify  the  nxn  matrix  X  =  [x2,  ....  xj,  Xj  €  Rn  and 
the  n2  vector  $  =  column  [xx,  .  .  .  ,  xn}.  If  we  then  identify  A^  with  -STX  (i.e., 
Aj  ■  -  I®ST)  and  A2<I>  with  -XS  (i.e.  A2  «  -S®I),  we  can  immediately  see  the 


equivalence  between  (4.11)  and  (4.12). 
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In  the  preceeding  sections  we  have  presented  an  algorithm  which  offers  some 
definite  advantages  in  computing  directly  the  feedback  gains  K  for  high  dimensional 
LQR  problems  such  as  those  arising  in  approximating  partial  or  delay  differential 
equation  control  problems.  As  we  shall  see  with  several  numerical  examples  in  this 
section,  it  can  substantially  outperform  standard  eigenvector  methods  on  such 
problems.  As  we  have  pointed  out,  a  fundamental  algebraic  operation  (in  both  the 
Chandrasekhar  update  (2.6),  (2.7)  and  in  the  reformulated  Smith  methods  (4.10b)) 
involves  computation  of 

(5.1)  L(I  -  r(A  -  BK))-1 

where  L  and  K  are  pxn  and  m*n  matrices  respectively.  Our  algorithm  uses  the 
Sherman-Morrison-Woodbury  formula  which  can  provide  significant  computational 
savings  when  m  and  p  are  small  compared  to  n.  For  systems  involving  sparse 
matrices  A  (a  frequent  occurence  in  many  approximation  schemes),  the  needed 
calculations  can  be  carried  out  quite  efficiently. 

We  further  note  that  the  Chandrasekhar  and  Newton-Kleinman-Smith 

components  as  formulated  in  our  algorithm  lead  to  ready  estimates  between  the  true 
gain  K  and  the  iterates  K;  in  terms  of  equation  errors  in  the  steps  being  performed. 

One  component  (the  variable  step  Smith)  of  the  algorithm  is  most  effectively 
carried  out  if  one  possesses  some  a  priori  knowledge  of  bounds  on  the  closed  loop 
eigenvalues.  If  the  closed  loop  eigenvalues  lie  close  to  the  imaginary  axis,  then 

convergence  in  the  Smith  method  can  be  very  slow.  Eigen  or  Schur  vector  methods 


[27],  [30]  arc  less  sensitive  in  this  regard.  For  low  order  systems,  the  Schur  vector 
approach  is  more  reliable  and  less  expensive  computationally  than  our  algorithm.  Our 
hybrid  algorithm  depends  critically  on  a  number  of  choices  (eg.,  stopping  criteria  in 
the  Newton-Kleinman  and  Smith  components,  stepsize  sequence  {r j)  and  iteration  count 


-19- 


sequence  {kj}  in  the  variable  step  component)  to  be  made  by  the  user  and  the  "best" 
choices  are  heavily  problem  dependent.  Hence  one  can  expect  our  hybrid  algorithm  to 
require  more  experimentation  and  fine  tuning  than  other  more  standard  methods. 

However,  as  we  shall  demonstrate  with  examples,  for  the  case  where  n  is  large 

compared  to  m  and  p,  it  can  offer  considerable  computational  savings  with  no  loss  in 
accuracy  over  the  methods  mentioned  above. 

We  have  tested  (and  are  continuing  our  efforts  in  this  direction)  our  hybrid 
algorithm  on  several  numerical  examples.  We  shall  report  on  just  two  of  these  here  to 
illustrate  our  findings.  All  our  computations  were  carried  out  in  double  precision  on 
an  IBM  3081  at  Brown  University.  We  gratefully  acknowledge  the  assistance  of  Yun 

Wang  in  our  carrying  out  of  the  extensive  computational  studies  reported  for  the 

boundary  flux  control  in  the  diffusion  equation  problem  of  Example  2  below. 

Example  1:  As  one  of  our  examples,  we  considered  an  example  (Example  6  of  [27]) 
which  Laub  used  to  test  his  Schur  based  methods.  The  system  is  the  n-dimensional 
system  of  (1.1)  with 
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which  leads  to  an  ill  conditioned  Riccati  equation.  This  problem  corresponds  to  one 
in  which  n  integrators  are  connected  in  series  with  a  feedback  controller  to  be 
applied  to  the  nth  integrator  in  order  to  stabilize  the  system.  Only  deviations  of  Xj 
from  the  origin  are  penalized  in  the  cost  functional.  The  true  optimal  gain  is  an  n 

vector  K  -  (K1 . Kn)  and  for  this  example  one  can  argue  that  K1  -  1.  In  (27), 

Laub  used  his  Schur  techniques  to  study  this  example  and  reported  difficulties  with 
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loss  of  accuracy  at  a  relatively  low  value  of  n,  n  ■  21.  We  carried  out  runs  with 
our  hybrid  algorithm  and  obtained  quite  favorable  performance.  Some  of  our  findings 
included: 

(a)  For  n  =  40,  we  used  the  Chandrasekhar  component  to  integrate  to  t{  =  100  and 
produce  an  initial  estimate  K.J  «  .99041,  which  when  used  in  the  Newton-Kleinman 
(fixed  step  size  r  =  .5  in  the  Smith)  produces  the  estimate  Kg  *  1.0  -  in  a  total  of 
2.93  seconds  of  CPU  time.  When  we  used  a  cruder  solution  in  the  Chandrasekhar 
component  (tf  ■  200  but  with  step  size  twice  that  in  the  first  run)  to  produce  K*  * 
.9394,  followed  by  the  N  -  K  (r2  =  .5,  r2  =  1.0  in  the  variable  step  Smith)  we 
obtained  Kg  =  1.0  -,  all  in  2.39  seconds. 

(b)  For  n  =  50,  we  produced  Kq  =  .9224745  at  tf  =  220  and  after  the  N-K-Smith 
(fixed  step  r  =  .5  in  the  Smith)  obtained  Kg  =  1.0000000003820  in  a  total  of  4.44 
CPU  seconds.  For  the  same  runs  with  variable  step  (rj  =  .5,  r2  ■  .7)  Smith  we 
obtained  a  Kg  as  above  with  3820  replaced  by  3817  in  a  total  of  4.31  seconds. 

(c)  We  compared  runs  with  the  Chandrasekhar  component  only  against  the  Potter 
method  for  n  *  10,  21,  40.  Obtaining  essentially  the  same  estimates  for  n  =  10  and 
21  (at  n  -  40,  the  Potter  degenerates  numerically  to  produce  useless  estimates)  we  had 
CPU  times  of  CHn=10  «  .753  seconds,  POTTn„10  «  .188  seconds,  CHn=21  -  1.52 
seconds,  POTTn=21  =  1.22  seconds,  CHn=40  =  4.35  seconds,  POTTn=40  -  6.81  seconds. 

We  found  for  this  example  that  the  eigenvector  methods  are  best  for  small 
n,  but  as  n  grows,  the  Chandrasekhar  alone,  and,  even  more  so,  the  hybrid  method 
will  out  perform  the  eigen-Schur  methods  in  both  accuracy  and  CPU  times.  A  more 
striking  demonstration  of  this  behavior  will  be  given  in  the  next  example. 

Example  2:  We  consider  the  linear  quadratic  regulator  problem:  minimize  the  cost 


functional 
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(5-2)  J(u)  =  J*(  |  Cz(t)  | 2  +  |u(t)  j2)dt 

subject  to  the  partial  differential  equation 

(5.3)  f-z(t,x)  “  —  z(t,x)  ,  x  €  (0,1) 

at  9x2 

z(0,x)  =  $(x) 

with  boundary  conditions 

(5.4)  z(t,0)  =  u(t)  and  |-z(t,l)  -  0 

where  c(-)  is  square  integrable  on  [0,1]  and 

Cz(t)  =  f1c(x)  z(t,x)  dx 
Jo 

We  can  discretize  or  approximate  (5.3)-(5.4)  using  the  standard  Galerkin  method  [2]; 
i.e.  the  approximating  solution  zN(t,x)  to  (5.3)-(5.4)  is  given  by 

(5.5)  zN(t,x)  =  £  W;(t)l  i(x)  ,  w.(t)  £  R1, 

i=0 

where  is  the  first  order  spline  defined  by 


Iftx)  = 


N(x  -  ^-), 

N 


n (i±L  -  x), 
N 


liill  < 


N 


x  <  L 
N 


L  <  x  <  iilii 

N  N 


otherwise 


and  zN(t,x)  satisfies 

(5.6)  J  |_  zN(t,x)<i>N(x)dx  -  -J  zN|_«/>N  dx  -  u(t)<Ji*tO) 

for  all  €  ZN=span  {£  . <j  )  . 

Then,  (5.6)  leads  to  the  nth  order  (n  -  N  +  1)  ordinary  differential  equation  for 
wN  -  col(w0 . wN)  ; 

(5.7)  QNwN(t)  -  -HNwN(t)  -  BNu(t) 


where 


and 


QN  =  i 


N 


H 


N 


N 


B 


N 


1/3  1/6  0 
1/6  2/3  1/6 
0 


1/6  2/3  1/6 

0  1/6  1/3 


with  QN  “ J1  Vjdx  . 


1  -1  0 

1  2  -1 
0 


-1  2  -1 
0  *1  1 


col  (1  0  ...  0)  . 


wi,h  HS  •  Hk '  hi  • 


For  computational  convenience,  we  change  coordinates  (for  fixed  N)  in  the 
system  (5.7)  by  x  =  QNwN  to  obtain  the  approximate  system 


x  =  -Hn(Qn)*1x 


Bn 


Thus,  in  (1.1)  we  have  A  -  -HN(QNyx,  B  -  -BN  and  C  -  CN(QNrx  where  CN  is  the 
vector  with  components  c^  -  c(x)lj(x)dx,  0  <  i  <  N. 

For  the  problem  in  this  example,  the  approximating  optimal  feedback 
operator  KN  is  given  [2]  by: 

Knz  -  Jj  kN(x)z(x)dx 

where  kN(x)  «  £  kjl.(x)  and  K  -  (k0,  ....  kN)  is  the  optimal  feedback  solution  in 

i=l 

the  problem  for  (1.1)  with  A,  B,  C  chosen  as  indicated  above.  We  note  in  this  case 
that  for  any  N  *  1,  A  has  only  one  unstable  eigenvalue  (zero),  (A,B)  is  stabilizable, 
and  (A,C)  is  detectable. 

For  the  special  case  when  c(x)  -  1,  we  find  C  -  (1,  ....  1)  <  r1x(n+1)  and 


hence  CA  -  0.  It  is  thus  easy  to  see  that  the  desired  solution  (K(t),L(t))  to  the 
Chandrasekhar  system  (1.3)  is  given  by 


where  k,i  are  scalar  functions  satisfying 


k  =  f  ,  k(0)  =  0 


i  «  -ik  ,  1(0)  -  1 

Therefore  we  find  kk  +  ii  =  0  so  that  kJ(t)  +  l2(t)  =  1.  We  thus  find  in  this  case 
that  k(t)  -*  1  as  t  -*  -®  and  hence  K  -  limK(t)  -  C.  For  this  case,  the 
Chandrasekhar  system  for  the  infinite  dimensional  LQR  problem  (5.2)-(5.4)  can  also 
be  analyzed  [17],  [36]  and  exactly  the  same  argument  as  above  shows  the  optimal 
feedback  gain  operator  is  given  by 
Kz  =  Jj  1  •  z(x)dx  . 

These  analytic  solutions  can  be  used  to  test  software  packages  and  approximation 
schemes  before  more  interesting,  analytically  intractable  examples  are  considered. 
Remark:  The  form  (5.7)  of  system  equations  appears  frequently  in  applications.  Thus 
the  critical  computational  factor  (5.1)  can  be  modified  so  that  one  can  avoid 
computing  A.  For  example,  in  this  case  it  has  the  form 
L(I  -  r  (-HQ’1  -  BK))'1 

(5.8) 

=  LQ(Q  +  rH  +  rBKQ)'1 

where  Q  +  rH  is  a  symmetric,  tridiagonal,  positive  matrix.  Thus  one  can  readily  use 
the  Cholesky  decomposition  algorithm  for  computing  LQ(Q  +  rH)'1  and  combine  this 
with  the  Sherman-Morrison-Woodbury  formula  (see  Remark  3  of  Section  2)  to 
efficiently  compute  the  critical  expression  (5.8). 

We  carried  out  extensive  computations  for  this  example  with  c(x)  -  1  +  x. 
We  compared  our  hybrid  method  to  the  Potter  algorithm  and  to  use  of  the 
Chandrasekhar  system  alone.  We  have  not  used  the  Laub-Schur  method  on  this 
example  since  we  felt  comparison  with  a  readily  available  (to  us)  Potter  package 
would  give  as  a  feel  for  the  relative  advantages  and  disadvantages  of  our  scheme 
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compared  to  eigen  and  Schur  vector  based  techniques.  (Analysis  and  computational 
experience  indicate  that  the  Potter  method  and  the  Laub-Schur  method  are  both 
0(NS)  with  the  latter  method  about  twice  as  fast  as  the  Potter  method.)  We  required, 
whenever  feasible,  the  same  level  of  accuracy  in  computation  of  feedback  gains  and 
compared  relative  CPU  times. 

In  studying  our  hybrid  scheme,  we  tested  numerous  sets  of  Smith 
acceleration  parameters  {rp,  {kp,  stopping  times  tf  in  the  Chandrasekhar  component 
and  error  stopping  criteria  in  both  the  Chandrasekhar  and  Newton-Kleinman-Smith 
components.  We  summarize  some  of  our  findings  to  date. 

In  Table  5.1  we  present  comparative  CPU  times  for  the  hybrid  scheme  vs. 
Potter  as  we  increase  N.  Recall  the  corresponding  finite  dimensional  approximation 
scheme  has  system  with  dimension  n  =  N+l.  In  all  of  the  runs  reported  in  Table  5.1, 
the  feedback  gains  for  the  hybrid  and  Potter  calculations  agreed  to  9  decimal  places 


t£ 

Hvbrid  (CPU  Sec.! 

Potter 

10 

.17 

.14 

20 

.31 

.81 

30 

.56 

2.45 

40 

.74 

5.49 

50 

.91 

10.71 

60 

1.09 

18.09 

70 

1.26 

27.97 

80 

1.43 

41.56 

100 

1.76 

120 

2.10 

140 

2.45 

160 

2.80 

Tafri?  5,1 

so  both  scemes  provided  accurate  solutions.  In  these  runs,  the  hybrid  scheme 
calculations  used  tf  =  2.2  (corresponding  to  h  -  .1)  with  |  L(-tf)  j  ■  10'3  in  the 
Chandrasekhar  component.  The  Newton-Kleinman  component  converged  after  4 
iterations  (i.e.  at  K4)  and  we  used  acceleration  steps  rx  =  1,  r2  -  10‘1,  rs  -  10'3, 


10 


-5 


Each  Smith  iteration  was  allowed  a  maximum  of  kj  -  50  per  value  of  rj 


although  in  most  cases  the  iteration  satisfied  a  convergence  criterion  before  this 
maximum  was  attained.  Careful  consideration  of  Table  5.1  reveals  that  the  hybrid 
scheme  is  clearly  O(N)  while  the  Potter  is  0(N3);  both  rates  are  to  be  expected  from 
our  earlier  observations  about  the  methods.  Note  that  at  N  *  80  the  hybrid  scheme 
is  more  than  25  times  faster  than  the  Potter  scheme  (with  comparable  accuracy,  of 
course). 

We  also  ran  the  hybrid  scheme  with  N  =  80  and  a  number  of  different 
fixed  acceleration  values  r  in  the  Smith  component.  The  same  Chandrasekhar 
component  prameters  as  reported  above  were  used.  Table  5.2  contains  relative  CPU 
times  as  well  as  an  indication  of  the  N-K  iterate  for  which  convergence  was 
achieved. 

In  Table  5.3  we  list  some  CPU  times  when  different  sets  of  acceleration 
parameters  {r^}  were  used.  Again  these  runs  were  for  N  =  80  with  the  same 


r 

Cm  (Sec.) 

Converged  N-K 

Gain 

5 

5.10 

K6 

1 

6.52 

K. 

10’1 

6.27 

K. 

10‘2 

7.40 

K. 

io-3 

10.10 

k. 

io-4 

12.06 

kb 

10'5 

10.07 

k4 

Table  5.2 


L  CPU  (Sec.) 

(10'\  10‘2)  6.25 

(1,  10-\  10’2)  4.88 

(1,  10*1,  10*2,  10’3,  10*\  10-5)  1.98 

(1,  10*\  lO"3,  lO-5,  lO'6,  lO’7)  1.61 


■Table  .5,3 

Chandrasekhar  solution  as  above.  All  of  the  converged  Newton-Kleinman  iterates  were 
after  6  steps  (i.e.  K6). 

Finally,  we  made  runs  (for  N  -  80)  to  find  the  best  results  that  the 


Chandrasekhar  algorithm  alone  (i.e.  accurate  integration  until  K(t)  -  K,  L(t)  -  0) 
could  produce.  The  best  results  we  were  able  to  achieve  yielded  an  accurate  value  of 
K  for  K(-tf)  with  tf  =  3.22  with  J  L(tf)  j  -  10*6  obtained  in  5.85  CPU  seconds. 

Based  on  our  computational  findings  for  the  above  two  examples  and  our 
experience  with  several  other  examples  for  infinite  dimensional  systems  (e.g.,  beams 
with  tip  bodies,  etc.),  we  are  quite  confident  that  the  hybrid  scheme  we  propose  in 
this  paper  can  be  profitably  used  with  a  number  of  large  scale  LQR  problems.  We 
are  currently  developing  a  rather  general  software  package  that  implements  the 
hybrid  scheme  in  a  manner  so  that  a  broad  range  of  problems  can  be  treated  in  the 
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